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ABSTRACT 


^This  document  presents  a  computational  system 
designed  to  remove  known  components  from  spectroscopic 
data.  The  system  realization  is  a  FORTRAN-77  program;  a 
listing  of  the  program's  source  code  is  included  in  the 
Appendix. 

The  program  first  performs  a  nine-point  quadratic 
least-squares  smoothing  function.  It  then  removes  the 
known  component  with  a  non-recursive  filter  based  on  a 
least-squares  prediction  of  the  observed  spectrum  and  a 
least-squares  prediction  of  a  computer  generated  synthe¬ 
tic  spectrum  of  the  known  component. 

The  program  was  tested  using  an  artificial  spec¬ 
trum  created  for  the  sole  purpose  of  testing  the 
program's  function. 

The  program's  function  was  verified  with  three 
segments  of  the  12  micron  band  of  ethane  taken  from  data 
recorded  on  the  Fourier  transform  spectrometer  at  the 
McMath  solar  telescope  at  the  Kitt  Peak  National  Observa¬ 
tory.  The  segments  are  located  near  804.0,  815.0,  and 
810.  s/cm'*^.  Plots  of  the  spectral  data  and  of  tlie  resul¬ 
tant  data  are  included. 
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CHAPTER  I 


INTRODUCTION 
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The  purpose  of  this  study  is  to  devise  a  technique 
whereby  any  known  component  of  a  given  spectrum  can  be 
removed  to  enable  study  of  the  underlyinq  spectral  data 
without  interference  from  the  unwanted  component.  The 
program  embodying  the  developed  technique,  called  FILTER, 
first  carries  out  a  quadratic  least-squares  smoothing 
operation  on  the  spectrum  being  studied.  Following  a  mild 
smoothinq  operation,  a  non- recu rsive  least-squares  filter 
is  applied  to  the  observed  spectrum.  This  filter  is 
constructed  from  the  predicted  synthetic  spectrum  of  the 
component  to  be  removed.  The  output  of  the  filter  may  be 
thought  of  as  an  "error  spectrum"  which  reflects  the 
least-squares  difference  between  the  observed  complete 
spectrum  and  the  synthesized  component  upon  which  the 
filter  is  based. 

Often  in  the  study  of  spectroscopic  data,  one  is 
interested  in  features  that  are  so  near  other  larger 
features  that  information  about  the  feature  in  question 
is  very  difficult  to  obtain.  Without  some  method  of 
removing  the  largo  features,  the  small  ones  arc  mostly 
inaccessible.  Of  course,  it  is  not  necessary  that  the 
desired  feature  be  a  small  one.  One  may  wish  to  separate 
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two  equivalent  sized  components  or  remove  a  known  small 
component  that  is  distorting  the  data  to  be  studied. 

For  instance,  in  the  study  of  ion  spectra,  any 
given  transition  is  extremely  weak  (  Sa  y  k  ally ,  1  9  8  4  ) . 
Unless  the  tra nsition  is  fortuitously  isolated,  it  may 
not  be  observable.  In  the  case  of  infrared  molecular 
spectroscopy  of  isotopic  molecules,  some  transitions  are 
separated  from  their  relatives  by  less  than  a  milli- 
wavenumber  (10~^  cm~M.  Clearly,  any  analysis  of  the 
spectrum  would  be  complicated  by  the  nearby  transitions. 
The  structure  of  planetary  atmospheres  can  be  determined 
through  the  use  of  spectroscopy,  but  many  tra nsitions- - 
and  consequently  many  atmospheric  const! tuents- -can  not 
be  identified  because  of  their  proximity  to  other  spec¬ 
tral  linos.  By  using  FILTER,  one  can  remove  selective 
portions  of  an  observed  spectrum  to  isolate  the  desired 
features  for  more  detailed  study. 

Another  potential  use  of  this  program,  and  the  one 
which  this  works  addresses,  is  in  the  study  of  hot-band 
spectra  which  normally  are  intermixed  with  the  cold-band 
spectra  and  consequently  are  difficult  to  pick  out  of  the 
data.  With  FILTER,  one  can  climinat'^  the  known  qround- 
state  transitions  from  the  obru^rved  spectrum  leaving  the 
smaller  hot-band  features  in  a  prominent  position  to  be 


studied . 


The  program  itself  is  a  simple,  yet  elegant, 
approach  to  a  complex  problem.  It  uses  a  nine-point 
least-squares  smoothing  routine  to  reduce  the  effect  of 
noise.  After  smoothing,  it  uses  the  calculated  synthetic 
spectrum  to  develop  a  predicted  component  of  the  observed 
spectrum  and  perform  a  non- rocu  r  siv  e  filtering  function. 
Once  the  smoothing/filtering  function  is  complete,  a 
simple  point-by-point  subtraction  leaves  the  error 
spectrum  in  the  resultant  data  file. 

The  data  files  must  first  be  processed  so  that 
they  represent  the  same  arrangement  of  inf ormation--data- 
point  density,  line  widths  and  so  forth  must  be  equiva¬ 
lent  or  the  results  will  be  meaningless.  Though  the 
manipulation  of  the  data  fil('s  required  to  achieve  such  a 
suitable  arrangement  for  inr>ut  to  the  program  represents 
an  additional  effort  in  the  analysis  of  the  spectrum,  the 
additional  information  available  in  the  residual  spectrum 
more  than  compensates  for  the  extra  work. 

The  discus  sion  in  this  work  is  geared  toward 
infrared  spectroscopy  and  the  tost  runs  were  accomplished 
using  infrared  spectra — specifically  of  the  band  of 
ethane,  but  the  program  is  c'asily  ada[:>table  to  any  spec¬ 
troscopic  data. 


CHAPTER  II 


BACKGROUND 

In  the  course  of  any  spectroscopic  analysis,  one 
will  encounter  several  concepts;  in  order  to  understand 
the  workings  of  FILTER,  one  needs  to  be  familiar  with 
some  basic  ideas. 

In  infrared  molecular  spectroscopic  theory  (Herz- 
berg,  194  5)  line  positions  and  intensities  are  the  pri¬ 
mary  data  used  to  obtain  physically  useful  and  meaningful 
results.  Once  position  and  intensity  have  been  accu¬ 
rately  established,  the  contribution  of  line  shape  to  the 
spectral  data  becomes  important  too.  In  addition  to  in¬ 
frared  theory,  one  should  understand  the  process  of 
deconvolution  (Blass  and  Halsey,  1981)  and  have  a  good 
gras[)  of  the  data  processing  techniques  used  to  manipu¬ 
late  the  spectra  (Pennington,  1  970).  The  use  of  the 

filtering  scheme  in  the  program  is  constrained  by  the 
digital  data-processinq  tecliniques  used  in  the  proaram 
(Bozic,  1  979).  When  one  is  knowledgeable  of  those  con¬ 
cepts,  the  validity  of  the  resultant  data  follows. 

Molecular  Vibration-Rotation  Theory 

The  physical  model  u.'-.ed  in  characterizing  molecu¬ 
lar  vribration-rotation  energy  levels  is  straightforward 
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but  complex.  Since  this  project  was  only  peripherally 
involved  with  the  concepts,  a  brief  summary  is  presented. 


One  considers  a  molecule  as  a  polyatomic  system  of 
N  nuclei  and  n  electrons.  The  system  then  has  3(N+n) 
degrees  of  freedom.  The  Born-Oppenheimer  approximation 
allows  one  to  separate  the  contribution  of  the  electronic 
motion,  so  only  3N  degrees  of  freedom  need  be  considered. 
Furthermore,  it  is  always  r'ossible  to  separate  the  three 
decrees  of  freedom  that  correspond  to  the  motion  of  the 
center-of-mass.  The  remaining  3N-3  degrees  of  freedom 
cannot  be  separated  easily,  but  by  assuming  the  molecule 
is  a  rigid  rotor,  one  can  assign  a  body-fixed  frame  of 
reference  and  separate  the  rotational  motion  from  the 
vibrational  motion  (Herzberg,  1  94  5). 


Ij  ine  Position 

This  body-fixed  frame  of  reference  allows  one  to 
assume  a  linear  combination  of  the  vibrational  and  rota¬ 
tional  Hamiltonians.  The  determination  of  molecular 
f'nerqy,  then,  is  a  combination  of  two  problems. 

The  solution  of  the  rotational  Hamiltonian  leads 
to  an  ox[)ression  for  the  rotational  energies  of  the 
mo  1  o'-u  1  e  : 


F  ( J,K)  =  B.7(J+1  )  +  (A-B)K 


when’  J  and  K  are  the  rotational  quantum  numbers  and  the 


rotational  constants  are 


D  -  - - - ,  -  - TT - .  V  ^  / 

Sn'^clg  Sn-^cl^ 

Similarly,  the  solution  of  the  vibrational  Hamiltonian 
leads  to  an  expression  for  the  vibrational  energies  of 
the  molecule: 


G(vj,V2, 


uij(v2  +  1/2)  + 

u>2(v2  +  1/2)  + 


(3) 


where  v^  are  the  vibrational  quantum  numbers  of  the  nth 
normal  mode  of  vibration  and  are  the  frequencies  of 

the  normal  vibrations  (Herzberg,  194  5). 

Unfortunately,  the  rigid-body  approximation  does 
not  give  sufficient  accuracy  for  high  resolution  spectro¬ 
scopy.  To  correct  for  the  non-rigidity  of  the  molecule, 
two  changes  are  needed.  First,  the  coefficients  are 
changed  slightly  from  the  equilibrium  values  used  in  the 
rigid-body  approximation  to  a  value  that  accounts  for  the 
average  motion.  Second,  some  additional  terms  are  added 
to  the  expressions  to  account  for  factors  such  as  anhar- 
monicity,  distortions  of  the  molecule  and  coriolis 
effects.  The  energy  expressions  are  rewritten: 

f''(J,K)  =  d'^J(J+1)  +  (A''-b'')k2  -  (At^.)Kl^ 
-dYj2(J+1)^  -  D,lij(J+l)K2  - 


(4) 


G(Vj,V2,  .  .  .)  =  +  1/2)  + 

i  k>i  l/2){v,^  +  1/2)  + 


where  the  v  superscript  indicates  an  average  value,  the 
D's  are  centrifugal  distortion  constants,  the  (A?^)  is  a 
Coriolis  constant  and  the  are  corrections  for  anhar- 
monicity.  The  total  vibration-rotation  energy  of  a 
molecule  is  then 


T  =  G(vj,V2,  .  .  .  )  +  f''(J,K) 


and  the  line  position  or  frequency  follows  from  the 
energy  difference  (Herzberg,  1945). 

The  literature  shows  that  a  good  study  of  most 
spectra  will  yield  accuracy  of  the  transition  coeffi¬ 
cients  to  within  approximately  .01%,  so  this  first  ele¬ 
ment  of  the  computer  program  has  a  good  foundation  and 
one  should  reasonably  expect  FILTER  to  give  accurate 
results  in  line  positions. 


liiH®  Intensity 

Line  intensity  is  not  so  simple  to  calculate  as 
line  position.  The  approximations  used  in  calculating 
line  intensity  are  more  susceptible  to  inaccuracy  than 
those  used  in  calculating  the  line  position.  These 
approximations  can  lead  to  proportionately  larger  errors 
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in  the  calculated  theoretical  intensity  than  should  be 
encountered  in  calculating  line  position  (Demtroder, 
1981  )  . 

Quantum  mechanically,  the  strength  of  a  given 
transition  is  proportional  to  the  square  of  the  transi¬ 
tion  probability  (Demtroder,  1981).  Unfortunately,  the 
literature  is  generally  unclear  in  describing  calculated 
line  strengths,  leaving  guess  work  or  t ria  1-and-error  as 
the  primary  tool  in  determining  what  has  been  included  in 
the  calculation  of  the  line  strength.  In  general,  the 
calculated  line  strength  must  be  scaled  by  some  conver¬ 
sion  factor  to  arrive  at  the  absolute  line  strength. 

Given  the  absolute  line  strength,  one  convolves 
the  calculated  line  strength  with  the  expected  line 
shape,  which  is  described  below,  to  get  the  absorption 
coefficient . 

The  calculated  absorption  coefficient  can  then 
used  in  the  expression  for  transmitted  intensity,  T: 

T  =  (7) 

where  is  the  incident  intensity;  a,  the  absorption 
coefficient;  P,  the  pressure  in  atmospheres;  and  1,  the 
f)ath  length  in  cm.  This  expression  for  T  is  derived  from 
Beer's  law  for  linear  absorption  (Demtroder,  1981).  The 
absorptance  obtained  is  then 


A 


1  .0  -  T 


1.0 


(8) 


- 


-aPl 


The  absorptance  space  obtained  with  this  expression  is 
the  most  accurate  space  available  for  analyzing  strongly 
absorbing  lines  in  a  spectrum. 

The  absorption  spectrum  is  sometimes  analyzed  in 
the  absorption  coefficient  space  obtained  by  a  Taylor 
series  of  the  expression  for  the  absorptance  which  yields 


(oPl  )2 


A  =  oPl  - 


In  the  first  approximation  of  this  Taylor  series,  when 
aPl  reaches  25%,  the  error  in  calculated  line  intensity 
is  at  least  3%.  Clearly,  whether  to  analyze  the  spectrum 
in  absorptance  space  or  in  absorption  coefficient  space 
is  dependent  upon  the  maximum  absorption  expected  in  the 
spectrum . 

Finally,  the  line  shape  has  not  yet  been  included 
in  the  calculation  of  line  intensity.  As  will  be  seen 
below,  the  line  shape  used  in  determining  the  theoretical 
spectrum  will  also  affect  the  intensity  of  the  peaks  in 
the  calculated  spectrum. 


Line  Shape 


The  expected  line  shape  is  equally  as  important  as 
the  line  strength  in  calculating  the  theoretical  spec¬ 
trum.  It  affects  the  calculated  line  intensity  during 
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the  convolution  process.  Furthermore,  the  shape  of  the 
expected  instrument  function  will  also  affect  the  calcu¬ 
lated  intensity. 

In  general,  the  intensity  of  a  given  line  will 
decrease  during  convolution  with  a  normalized  profile 
because  the  convolution  increases  the  width  of  the  line. 
The  overall  problem  of  selecting  a  line  shape  and  con¬ 
volving  that  line  shape  with  the  stick  spectrum  contains 
many  factors  which  will  alter  the  intensity  of  the  theo¬ 
retical  spectrum. 

The  actual  line  shape  is  taken  to  be  a  combination 
of  several  factors.  The  natural  line  profile  is  a  Lorent- 
zian.  It  is  on  the  order  of  a  few  kHz  and  for  the  pur¬ 
poses  of  this  study  it  can  be  approximated  by  a  delta 
function  (Demtroder,  1981).  The  doppler  line  width  is 
represented  by  a  Gaussian  whose  full  width  at  half  maxi¬ 
mum  (FWflM),  Ap,  is  a  function  of  the  temperature  of  the 
sample : 

Ap  -  7.16  X  lO""^  v^yT/M  (cm"^)  (10) 

where  T  is  the  temperature  in  °K  and  M  is  molecular  mass 
in  amu.  Likewise,  the  observed  line  [irofile  includes  a 
contribution  due  to  the  pressure  of  the  sample.  The 
pressure  broadening  can  be  represented  by  a  Lorentzian 
with  FWHM 


0.250P  (cm 


(11) 


where  P  is  the  pressure  in  atmospheres  (Blass,  1985). 

Both  the  doppler  width  and  the  pressure  width  must 
be  taken  into  account  when  calculating  the  expected  line 
shape  of  the  sample  to  be  studied.  The  resultant  line 
shape  is  then  convolved  with  the  calculated  stick  spec¬ 
trum  to  get  the  theoretical  absorption  coefficient  spec¬ 
trum.  In  absorptance  space,  the  calculated  spectrum  can 
be  expressed 


A  =  1  -  ® 


(12) 


where  a  indicates  a  convolution  function,  D  is  the  dop¬ 
pler  width  profile,  P  is  the  pressure  width  profile,  S  is 
the  absolute  line  strength,  and  PI  is  the  pressure-path 
length  product  of  the  sample.  The  absorption  coeffi¬ 
cient,  a,  then  is  the  convolution  of  D  and  P  times  S. 


Deconvolution 


In  addition  to  the  "real"  line-shape  factors  men¬ 
tioned  in  the  previous  section,  one  must  also  consider 
the  instrument  function  when  working  with  spectral  data. 
After  the  theoretical  absorptance  spectrum  has  been 
determined,  the  expected  instrument  function  is  convolved 
with  the  theoretical  absorptance  spectrum  to  generate  a 


calculated  synthetic  spectrum  containing  known  compo¬ 
nents  . 


The  observed  spectrum  necessarily  contains  the 
instrument  transfer  function,  which  broadens  the  observed 
line  width,  effectively  reducing  the  resolution.  Decon¬ 
volution  is  a  technique  that  can  be  used  to  increase 
apparent  resolution  or  reduce  required  observing  time  by 
removing  the  instrument  function  from  the  observed  data. 
The  accuracy  of  the  deconvolution  process  is  highly  de¬ 
pendent  on  a  signal  to  noise  ratio  of  100  or  better 
(Blass  and  Halsey,  1981). 

Fi Itering 

The  digital  spectral  data  may  be  filtered  either 
as  it  is  taken  or  it  may  be  recorded  in  raw  form  and 
filtered  in  the  process  of  analysis.  The  latter  approach 
offers  the  advantage  of  safety  and  flexibility.  If  the 
data  were  to  be  filtered  as  it  wore  taken,  it  could  be 
catastrophically  altered  as  a  by-product  of  the  fil¬ 
tering.  By  recording  the  data  before  filtering,  one  may 
always  start  over  if  non-physical  results  are  achieved. 
Furthermore,  more  than  one  filtering  method  may  be  used 
if  the  results  are  not  suitable. 

One  may  work  with  either  continuous-time  or 
discrete -time  data.  Continuous-time  refers  to  a  data  set 
in  which  the  independent  varial>le  may  assume  a  continuous 


range  of  values  while  discrete-time  implies  that  the 
independent  variable  is  restricted  to  a  finite  discrete 
set  of  values.  A  discrete-time  signal  is  often  referred 
to  as  a  sampled  data  signal  (Freeman,  1965). 

The  digital  signal  may  have  data  points  taken  at 
any  of  the  discrete  values  of  the  independent  variable. 
These  values  are  usually  equally  spaced,  however  equal 
spacing  is  not  a  requirement.  The  data  processing  is 
considerably  simpler  with  equal  spacing  (Bozic,  1979). 
In  the  discrete-time  system,  the  independent  variable 
time,  t,  is  taken  as 

t=kT 

where  T  is  the  difference  between  any  two  adjacent  data 
points  and  k  corresponds  to  the  kth  data  point.  The 
notation  used  in  this  work  is  x(k)  and  y(k)  for  the  kth 
value  of  the  independent  and  dependent  variables  respec¬ 
tively. 

Tlie  process  of  filtering  involves  predicting  a 
data  point  based  on  a  number  of  input  data  points.  It 
usually  consists  of  applyina  some  operation  such  as  delay 
or  multiplication  by  a  constant  to  the  input  data  points. 
If  the  filter  uses  past  output  as  input,  it  is  called 
recursive,  while  if  it  uses  no  past  output,  it  is  called 
non-recursive  or  moving  average.  Tlie  recursive  filter  has 


some  problems  with  stability  and  was  not  considered  for 
this  study  (Bozic,  1979). 


The  general  representation  of  a  digital  filtering 
scheme  is 


y(k)  +  bjy(k-l)  +  .  .  .  +  bj^y(k-m)  = 

aQx{k)  +  a]^x(k-l)  +  .  .  .  +  aj^x(k-n) 


(13) 


where  x(k)  represents  the  kth  input  point;  y(k),  the  kth 
output  point;  and  the  s  and  bj^'s  are  weighting  factors 
with  b^  taken  as  unity  (Boz ic,  1979).  This  representa¬ 
tion  means  that  at  time  k  (t=kT)  the  output  of  the  filter 
is  the  current  input  plus  some  linear  combination  of 
previous  inputs  and  outputs.  In  the  case  of  using  re¬ 
corded  data,  one  may  also  include  future  inputs  and 
outputs.  For  the  non-recursivc  filter  all  bj^,  m>0,  are 
zero.  Figure  lisa  generalized  depiction  of  a  recur¬ 
sive  filter  while  Figure  2  illustrates  a  non-recursive 
filter. 

Digital  Processing 


Use  of  poor  technique  in  digital  data  processing 
can  lead  to  catastrophic  alteration  of  the  data;  hence, 
one  must  be  alert  not  to  generate  excessive  errors  in 
processing  the  data.  There  are  two  factors  in  digital 
data  processing  that  can  contribute  to  the  error:  quan¬ 
tization  and  finite  word  length. 


^  J 


Figure  2.  Diagram  of  a  non-rccursivc  filter, 


Quantization  error  results  from  trying  to  describe 
infinit-^  precision  data  with  finite  precision.  The  posi¬ 
tion  o  1-  value  of  a  data  point  will  never  precisely  coin¬ 
cide  with  the  discrete  value  available  to  describe  it. 
Furthermore,  there  is  some  quantization  error  involved  in 
the  original  measurement  of  the  spectrum,  but  that  is 
beyond  the  scope  of  this  study.  The  slight  inaccuracies 
involved  in  quantization  error  will  appear  as  noise  in 
tliG  spectrum. 

Finite  word  length  errors  will  result  in  a  loss  of 
accure.  ey  and  generate'  aj'parent  random  noise  when  the 
result  of  any  operation  exceeds  the  precision  of  the 
com[)Ut"r.  As  a  sim[)le  example,  assume  the  computer  can 
r>tor<’  only  five  digits  plus  the  ox[)onent.  Ttic'n  let  x  = 
1  .4  34  d  X  10^  and  let  y  =  1.73  32  x  10^.  Wlien  x  is  added 
to  y,  tlie  result  is  1  8,766,^00,  hut  the'  I'^sult  will  be 
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stored  as  1.8766  x  10  and  the  information  that  the 
number  is  actually  larger  by  500  will  be  lost. 

These  errors  are  cumulative  and  the  effect  can  be 
much  larger  than  the  one  given  in  this  example.  Two 
procedures  are  useful  in  minimizing  the  effect  of  finite 
word  length  errors:  avoiding  subtraction  of  a  large 

number  from  another  large  number  and  multiplying  by  fac¬ 
tors  less  than  one  whenever  possible  (Pennington,  1970). 


1  7 


CHAPTER  III 


DESCRIPTION  OF  THE  PROGRAM 

FILTER  is  written  in  FORTRAN-77.  An  attempt  was 
made  not  to  use  machine-specific  code  so  that  the  program 
can  bo  used  on  other  machines  with  only  minor  alteration, 
however,  the  data  files  are  identified  to  the  program 
through  interactive  input/output.  If  interactive  input/ 
output  is  not  possible,  the  source  code  must  be  modified 
to  provide  identification  of  the  spectra  to  be  used.  The 
source  code  for  FILTER  is  listed  in  the  Appendix. 

The  spectrum  files  are  formatted  in  double  preci¬ 
sion  format  with  each  point  written  as  D20.12.  The 
starting  frequency  of  the  spectrum  is  the  first  entry  in 
the  file  and  the  step  size  is  the  second.  Following  the 
starting  frequency  and  the  stop  size,  the  spectral  inten¬ 
sities  are  v/ritten  to  the  file.  The  file  must  be  written 
for  random  access  input/output.  If  the  data-point  density 
and  starting  frequency  of  the  calculated  spectrum  and  the 
observed  spectrum  are  not  identical,  the  calculated  spec¬ 
trum  must  be  altered  in  order  to  coincide  with  the 
observed  spectrum. 

The  program  uses  a  buffc'rcd  input/output  capabil¬ 
ity  to  handle  large  data  files  without  using  a  prohibi¬ 
tive  amount  of  memory.  It  currently  works  with  a  100 
point  data  set,  buffered  by  10  points  on  either  end  to 
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rocluc'^  i!!\riul5;r->  effoct^^.  Nonnally,  tho  pronrnm  gots  the  10 
buffer  ncDints  from  the  spectrum  itself,  using  the  10 
points  ^dj,^cer;^  to  the  data  set  in  use.  If  the  data  set 
in  use  is  the  first  data  set  in  tho  spectrum.  Subroutine 
RBt'f'Fl  creates  a  buffer  of  10  [)oii'.  ts  b  v'  applying  an 
e  >;  [Ton  e-n  t  i  a  1  d'''cay  to  the  first  point  in  the  data  set, 
after  considering  tlie  slope  of  th.o  spectrum  at  that 
point.  Subroutine  RBUI''F2  [K'r  forms  tlie  same  f’.inction  on 
the  trailing  end  of  tho  last  data  sc-t  in  the  spectrum. 
It  should  be  noted  that  100  point  (bit  a  sots  ai'o  unrealis¬ 
tically  small  for  most  spectroscopic  a  ppl  i  cations.  The 
sise  VMS  cl'ioscn  initially  to  minimise  the  cemfuitor  time 
used  in  tho  dcw(dopm(?nt  phase';  tlio  exte^nsion  to  larcer 
blocks  is  straightforvviard. 

Subreutino  REUFF  roads  the  spectrum  into  memory 
using  the  120  point  data  set;'.--l()0  points  ol  working  data 
and  20  points  of  buffer  data.  To  illustrate  the  function 
of  the  program,  one  data  set  will  hi'  tracked  Ihrounh  its 
proee  s  s i ng . 
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less  than  unity.  If  necessary,  LSQ  calls  subroutine 
INVERT.  INVERT  pivots  the  matrix  to  reduce  error  propa¬ 
gation  by  insuring  that  all  multiplications  are  by  fac¬ 
tors  less  than  one. 

The  next  step  is  the  actual  filter  application. 
Subroutine  FILTER  uses  a  non-recursive  filter  generated 
from  a  least-squares  prediction  for  two  points  on  either 
side  of  the  nth  point.  The  filter  weights  the  predic¬ 
tions  for  the  nth  point  based  on  the  number  of  data 
points  separating  the  prediction  from  the  nth  point.  The 
predictions  two  points  away  are  weighted  by  a  factor  of 
0.1,  those  one  point  away  are  weighted  by  0.2  and  the 
prediction  at  the  nth  point  is  weighted  0.4.  Figure  3 
depicts  the  filtering  scheme  used  for  the  interior 
points . 

Next,  FILTER  shifts  one  point  to  the  right  and 
repeats  the  process  for  the  next  point,  producing  a 
"moving  average."  It  should  be  noted  that  the  end  points 
of  the  data  set  have  a  different  weighting  scheme  which 
is  enumerated  in  the  source  code,  liowever,  points  calcu¬ 
lated  with  that  different  weighting  scheme  are  all  dis¬ 
carded  when  the  buffer  points  are  stripped  from  the  data 
set. 


After  the  filter  application,  subroutine  WBUFF 
stri[)s  the  2  0  buffer  points  from  the  data  set  and  sequen¬ 
tially  writes  the  100  working  data  points  to  a  temporary 


holding  file.  Once  the  working  data  set  is  written  to 
its  holding  file,  its  existence  as  a  data  set  is  com¬ 
plete  . 

The  program  loops  back  to  subroutine  RBUFF  to 
begin  the  process  for  the  next  data  set  in  the  spectrum. 
The  cycle  continues  until  there  are  no  more  100  point 
data  sets  remaining  in  the  spectrum. 

When  processing  is  complete  for  the  observed  spec¬ 
trum,  the  program  begins  processing  the  calculated  spec¬ 
trum,  repeating  the  entire  process  that  was  performed  on 
the  observed  spectrum  for  the  calculated  spectrum. 

With  each  spectrum  in  its  respective  holding  file, 
subroutine  MINUS  reads  the  processed  observed  spectrum 
and  the  processed  calculated  spectrum  into  memory  and 
subtracts  the  calculated  spectrum  from  the  observed  spec¬ 
trum  point-by-point  to  got  the  residual  error  spectrum. 

Finally,  before  the  program  terminates,  an  option 
to  subtract  an  additional  calculated  spectrum  is  offered. 

A  cursory  error  trapping  function  is  included  in 
subroutine  ERROR.  If  an  error  is  detected  by  RBUFF  while 
reading  a  spectrum  file,  ERROR  issues  a  warning  that  the 
error  was  encountered.  It  should  be  noted  that  ERROR 
does  not  check  the  validity  of  the  input  data  and  the 
program  does  not  terminate  if  an  error  is  encountered. 
The  sole  purpose  of  ERROR  is  to  caution  the  user  that  his 
residual  spectrum  is  probably  inaccurate.  If  ERROR 


CHAPTER  IV 


RESULTS 

Testing  Using  Simulated  Data 

Since  the  concept  of  the  program  is  new,  testing 
presented  a  problem.  There  are  no  known  spectra  with 
which  to  compare  the  results  to  determine  that  a  valid 
solution  has  been  achieved.  Consequently,  a  "known" 
spectrum  was  created  for  testing  purposes. 

First  two  individual  stick  spectra,  called  DA  and 
DB,  were  created  with  a  pseudo-random  number  generator. 
Line  positions  were  randomly  determined;  then  a  random 
intensity  was  matched  against  each  line  position.  A 
third  spectrum,  D2,  was  created  by  combining  the  first 
two.  Finally,  each  spectrum  was  convolved  with  a 
Gaussian  whose  FWHM  was  seven  data  points  to  simulate  the 
combined  effects  of  the  natural  line  width,  Doppler 
broadening,  pressure  broadening,  and  the  instrument  func¬ 
tion.  The  convolved  spectra  are  called  GA,  GB  and  G2 
respectively . 

The  final  G2  spectrum  then  contained  two,  and  only 
two,  jsnown  components.  If  FILTER  was  applied  to  remove 
GA  from  G2,  the  resultant  should  have  been  GB.  The 
results  are  illustrated  in  Figures  4-7.  A  further  check 
was  made  by  removing  the  calculated  GB  from  the  residual 


Figure  6.  Plot  of  G2 


of  G2  -  GA  (for  convenience,  the  act  of  removing  one 
spectrum  from  another  will  be  denoted  by  The  final 

result  was  an  error  spectrum  whose  peak  intensity  was  on 
the  order  of  1.0  x  10”^. 

To  complicate  the  picture,  a  "noise"  spectrum  was 
generated,  added  to  the  stick  spectrum  D2  and  convolved 
with  the  Gaussian.  The  noise  spectrum  was  generated 
using  a  pseudo- random  number  generator  similar  to  the  one 
used  to  generate  DA  and  DB.  The  resultant  spectrum  was 
scaled  to  give  a  signal-to-noise  ratio  of  approximately 
20:1.  After  removing  GA  from  the  new  g2,  the  result  was 
GB  to  within  the  noise.  G2  -  GA  is  illustrated  in 
Figure  8. 

The  final  step  in  the  testing  phase  was  to  add  a 
single  lino  to  the  D2  stick  spectrum  that  contained  the 
additional  noise  spectrum.  The  intensity  of  the  line  was 
approximately  10%  of  the  average  local  intensity  and 
approximately  twice  the  average  local  noise  level.  This 
new  G2  spectrum  was  then  convolved  with  the  Gaussian,  and 
FILTER  was  used  to  remove  GA  and  GB.  The  resultant 
spectrum  faithfully  retained  the  added  lino  and  removed 
all  other  lines  above  the  noise.  See  Figure  9. 

Actual  Data  Runs 


After  validating  the  basic  concept,  FILTER  was 
used  to  remove  the  calculated  band  of  ethane  from  the 


observed  ethnne  spectrum  at  three  different  frequencies. 
The  frequencies  were  chosen  to  represent  three  different 
levels  of  activity:  one  with  low  activ'ity  near  815.0 

cm“^,  another  with  some  activity  hut  little  expected  hot- 
band  activity  near  804.0  cm“^,  and  a  final  one  with  a 

great  deal  of  cold-band  activity  and  hot-band  activity 
near  819.5  cm~^. 

The  segments  of  the  observed  ethane  spectrum  used 
were  recorded  on  the  one-meter  Fourier  transform  spectro¬ 
meter  (FTS)  at  the  McMath  solar  telescope  at  Kitt  Peak 
National  Observatory  (Jennings,  198  4).  This  is  similar 
to  data  taken  by  Daunt  et  aj^.  (  1  984)  except  the  resolu¬ 
tion  was  doubled  to  approximately  2.5  milli-wa venumbers 
by  modifying  the  FTS  to  double  pass  the  interferometer. 
The  calculated  stick  spectrum  used  is  taken  from  data  by 
Abakan  ejt  (1  983). 

In  computing  the  theoretical  spectrum,  it  was 
found  that  a  factor  of  12.1  x  10  cm  atm  was  neces¬ 
sary  to  convert  the  relative  intensities  given  in  the 
spectral  catalog  to  absolute  intensities.  The  doppler 
width  was  calculated  to  be  1.8  mil  1  i -wavenumbe rs  and  the 
pressure  width,  0.5  milli-wavenumbers.  The  instrument 
function  for  the  FTS  was  taken  to  be  a  sine  function  with 
FWllM  of  3.0  milli-wavenumbers.  The  o)5served,  calculated 
and  resultant  spectra  are  illustrated  in  Figures  10-18. 
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Calculated  v„  band  of  ethane  near  804.0  cm 


Portion  of  the  observed  spectrum  of  ethane  near  815.0  cm 


Calculated  band  of  ethane  near  819.5  cm 


As  an  additional  cross  check  of  the  validity  of 
the  resultant  data,  one  very  small  line  was  removed  from 
the  calculated  stick  spectrum  at  819.66089  cm~^.  The 
continuous  spectrum  calculated  from  this  new  stick  spec¬ 
trum  is  shown  in  Figure  19,  with  the  position  of  the 
removed  line  indicated  by  an  asterisk  {*).  The  error 
spectrum  produced  by  using  the  new  continuous  spectrum  in 
FILTEJR  is  identical  to  the  error  spectrum  produced  by 
using  the  original  stick  spectrum  except  for  the  reap¬ 
pearance  of  the  single  line.  The  line  is  identified  with 
an  asterisk  (*)  in  Figure  2  0. 


CHAPTER  V 


DISCUSSION 

The  spectra  used  in  the  testing  phase  are  fully 
known  since  they  were  created  specifically  for  testing. 
Consequently,  evaluation  of  the  program  function  is 
fairly  objective.  G2  with  GA  and  GB  removed  is  effec¬ 
tively  zero.  Furthermore,  when  G2  contains  features  not 
in  GA  or  GB,  those  extra  features  remain  when  GA  and  GB 
are  removed.  Hence,  one  can  say  the  procedures  employed 
by  FILTER  are  correct  for  the  intended  purpose. 

When  considering  real  spectral  data,  the  interpre¬ 
tation  of  the  results  is  definitely  more  subjective; 
there  is  no  reliable  predictor  of  the  residual  spectrum. 

Before  evaluating  the  program's  effectiveness  with 
observed  spectra,  one  must  ensure  that  the  observed  spec¬ 
trum  and  the  theoretical  spectrum  truly  represent  the 
same  data.  That  is,  each  data  point  in  the  observed 
spectrum  must  be  exactly  equivalent  to  the  corresponding 
data  point  in  the  calculated  spectrum.  Failure  to  pro¬ 
perly  match  the  spectra  will  generate  some  oscillatory 
behavior  or  large  negative  excursions  by  mismatching  peak 
positions.  Neither  condition  is  desirable  and  any  nega¬ 
tive  absorption  is  clearly  non-physical.  To  ensure  that 
the  data  is  in  fact  equivalent,  the  starting  frequencies 


must  exactly  coincide  and  the  data-point  density,  or  the 
step  size,  must  also  be  the  same. 

Furthermore,  when  converting  the  calculated  stick 
spectrum  to  a  continuous  theoretical  spectrum  as  in  Equa¬ 
tion  (12),  Page  11,  the  line  shape  chosen  to  convolve 
with  the  stick  spectrum  will  affect  the  peak  intensities 
of  the  calculated  spectrum.  A  normal! zed  line  shape 
preserves  the  total  energy,  but  in  doing  so  it  neces¬ 
sarily  reduces  the  peak  intensity  by  spreading  the  energy 
over  a  wider  frequency  range.  In  the  present  work,  the 
absorption  coefficient  was  calculated 


a  =  (0  a  P)Sq 

where  is  the  absolute  line  strengtli  determined  from 

data  by  Atakan  et  al_.  (1983)  as 

Sq  =  7.1  X  10"^  S,^ 

where  Sr  is  the  published  relative  line  strength. 

Using  data  supplied  by  Jennings  (1  984)  of  1.5  torr 
pressure,  296°K  temperature,  and  150  cm  path  length,  it 
was  found  that  the  intensities  of  the  calculated  spectrum 
were  too  low  by  a  fcTctor  of  1/3  to  effectively  remove 
their  contribution  from  the  observed  spectrum.  In  order 
to  bias  the  intensities  of  the  calculated  spectrum  to 
match  those  of  the  observed  s[iectrum,  a  factor  of  12.1  x 


10”'^  atm~'^  was  used  in  converting  the  calculated 

relative  intensities  to  absolute  intensities. 

The  discrepancy  between  the  expected  and  actual 
conversion  factor  can  be  ascribed  to  inaccuracy  in  the 
measurement  of  the  pressure  of  the  sample  (Jennings, 
1  984).  If  one  assumes  the  given  value  of  7.1  x  10”^  is 
correct,  the  discrepancy  corresponds  to  an  actual  pres¬ 
sure  of  2.6  torr.  Additionally,  the  discrepancy  may  be  a 
result  of  inaccuracies  in  establishing  the  base  line 
(Blass,  1  98  5).  Obviously,  any  coml:)ination  of  the  two 
errors  could  account  for  the  discrepancy. 

The  natural  question  arising  is  "IJow  does  one 
decide  that  the  calculated  spectrum  is  a  good  predictor 
of  the  Ijnown  part  of  the  observed  spectrum?" 

By  lodging  at  a  section  of  the  spectrum  with 
relatively  little  activity,  one  may  mal<e  an  educated 
guess  at  the  reliability  of  the  residual  spectrum.  Pre¬ 
sumably,  most  of  the  lines  will  be  completely  eliminated 
by  the  filter,  leaving  noise.  The  spectrum  near  815.0 
cm~^  demonstrates  that  FILTER  will  remove  those  compo¬ 
nents  when  there  is  not  much  underlying  activity. 

For  the  more  active  sections  of  the  spectrum,  a 
comparison  of  the  calculated  spectrum  with  the  observ'^ed 
spectrum  can  be  used  to  determine  tliat  a  given  lino  has 
been  completely  removed  from  the  residual  spectrum.  In 
such  a  case,  care  must  be  exercised  not  to  confuse  an 


underlying  element  of  the  observed  spectrum  that  should 
properly  exist  in  the  residual  spectrum  with  a  portion  of 
the  calculated  spectrum  that  was  not  completely  removed 
by  FILTER.  The  spectrum  near  804.0  cm~^  includes  a 
segment  of  fairly  high  spectral  activity  followed  by  a 
segment  of  relatively  little  activity  showing  that  FILTER 
will  remove  the  known  components  in  a  region  of  high 
spectral  activity  without  altering  the  underlying  data. 
The  spectrum  near  819.5  cm~^  is  an  extremely  busy  portion 
of  the  ethane  spectrum,  complete  with  hot-band  activity 
underlying  the  ground-state  transitions,  yet  FILTER 
effectively  removed  the  ground-state  transitions,  leaving 
what  is  presumed  to  be  hot-band  transitions.  Analysis  of 
the  residual  spectrum  is  left  for  further  study. 

Whether  to  deconvolve  the  data  before  applying 
this  program  is  a  question  to  be  answered  before  the  data 
is  taken.  Expected  resolution  of  the  experimental  data  or 
an  unreasonably  long  observation  are  factors  to  consider, 
but  the  use  of  FILTER  is  not  dependent  on  the  existence 
of  deconvolved  data.  However,  the  effect  that  FILTER  may 
have  on  the  deconvolution  process  has  not  yet  been  inves¬ 
tigated,  and  the  very  strong  sensitivity  of  the  deconvo¬ 
lution  process  to  the  signal-to-noise  ratio  of  the  data 
suggests  that  deconvolution  should  be  performed  before 
the  algorithm  in  FILTER  is  applied. 
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PROGRAM  LISTING 


. . . . . . . . . . . . . 

C*  PROGRAM  FILTER 

. . . . . . . . . . . . 

c 

C  THIS  PROGRAM  ACTS  AS  A  FILTER  TO  REMOVE  KNOWN  COMPONENTS  FROM  A 
c  DicitAi.  sPEcrr«jM.  IT  sMOoncs  III:  siicira  u;,ing  a  yoAORAiic 
C  LEAST  SQUARES  FIT  WHERE  NFlI«2tl  IS  THE  NUMUER  OF  POINTS  IN  THE  LEAST 
C  SQUARES  FIT,  THEN  IT  USES  THE  TWO  POINTS  ON  EITHER  SIDE  OF  A  POINT  AND 
C  TT€  SAME  POINTS  IN  THE  CALCULATED  SPECTRUM  TO  MAKE  A  mEDICTION  FOR  EACH 
C  POINT  IN  TTE  SPECTRUM.  TO  ACCOMODATE  LARGE  SPECTRA,  A  BUFFERING  ROUTINE 
C  IS  USED  TO  WORK  WITH  SECTIONS  OF  THE  Sf^CTRA  INDIVIDUALLY.  IT  CURRENTLY 
C  USES  100  DATA  POINTS  BUFFERED  BY  10  DATA  POINTS  ON  EITHER  END. 

C 

DIMENSION  XINT1(120),XINT2(120),A1(120),A2(120) 

DIMENSION  B1(120),B2(120),C1(120),C2(120),FILE1(6),FILE2(6) 

DOUBLE  PRECISION  XINTl,  XINT2,  A1 ,  A2,D1 ,02,01 ,C2 
COMMON  NDATA,NGIT,NBUFF,NFLAG 
C 

C  NDATA  IS  THE  DATA  SET  SIZE  PLUS  THE  TOTAL  NUMBER  OF  BUFFER  POINTS 
C 

NDATA=120 

C 

C  NBUFF  AND  TTLAG  ARE  FLAGS  USED  BY  SUBROUT  INES 
C 

NDUFP=0 
NFL AG =0 
C 

C  ISI2E  IS  THE  ESTIMATED  SIZE  OF  THE  SPECTRA/10  +1.  IT  IS  USED 
C  AS  A  LIMIT  TO  PREVENT  AN  INFINITE  LOOP 
C 

ISIZE=11 

C 

C  NFIT»2+1  =  THE  NUMBER  OF  POINTS  USED  IN  THE  SMOOTHING 
C 

NFIT=4 

WRITE  (5,100) 

100  FORMATC  input  FILESPEC  FOR  THE  MAIN  SPECTRUM') 

READ  (5,200)  FILEl 
200  FORMAT  (6A5) 

DO  20  I=1,ISIZE 

CALL  RBUFF(XINT1, FILEl, INI) 

IF  (INT+IO.LT. NDATA)  N0ATA=INT+10 
CALL  L5Q(XINT1,A1,B1,C1) 

CALL  FILTER(XINT1,A1,B1,C1) 
call  WBUFF(XINT1, FILEl, INT) 


IF  (NFLAG.EQ.l) 

GO 

TO  30 

20 

CONTINUE 

30 

WRITE  (5, AOO) 

AOO 

FORMAT  (•  input  FILESPEC 
READ  (5,200)  FILE2 

OF 

TIE  SPECTRUM  TO  DE  SUBTRACTED') 

c 

C  REINITIALIZE  THE  FLAGS 
C 

NDUFFrO 

NFL AG ^0 

DO  AO  I=1,ISIZE 

CALL  RBUFF(XINT2,FILE2,INT) 

IF  (INT+IO.LT. NDATA)  NDArA=INT4lO 
CALL  LS0(XINI2,A1 ,U1,C1) 

CALL  FILTER(XINT2,A1 ,B1,C1 ) 

CALL  WBUFF(XINT2,FILE2,INT) 

IF  (NFLAG.EQ.l)  GO  TO  50 
AO  CONTINUE 

50  CALL  M]NUS{F]LE1,FILE2,INI) 


53 


c 

C  OFFER  OPTION  TO  REMOVE  MORE  FROM  Tl€  SPECTRUM 
C 

WRITE  (5,500) 

500  FORMAT  ('  IS  THERE  ANOTHER  SPECTRUM  TO  BE  SUBTRACTED?  I=YES') 

READ  (5,600)  LOOP 
600  FORMAT  (II) 

IF  (LOOP.EQ.I)  GO  TO  30 
WRITE  (20,300)  XINTl 
300  FORMAT  (020.12) 

END 

»*♦**»****»»*»*»*»»»»»»»  •»•»»»•»»•••»«»»•»•*»•»•»**••*••«**•******• 

SUBROUTINE  LSQ(Y,A,a,C) 

0»***»*»»********»*»*******»»**»»»*»***»»*« ••**•*»••«*»* ********************** 

C 

C  THIS  SUBROUTINE  MARES  A  QUADRATIC  LEAST  SQUARES  FIT  TO 
C  A  DIGITAL  SPECTRUM.  NDATA  IS  THE  NUMBER  OF  DATA  POINTS. 

C  NFIT*2+1  IS  THE  NUMBER  OF  POINTS  IN  THE  LEAST  SQUARES  FIT 
C 

COMMON  ^DATA,NFIT 

DIMENSION  MATRIX(3,4),  Y(NDATA),  A(NDATA),  B(NDATA),  C(NDATA) 

DOUBLE  PRECISION  MATRIX,  Y,  A,  B,  C,  MULT 
DO  10  N=l, NDATA 
C 

C  INITIALIZE  THE  MATRIX  ELEMENTS 
C 

SXUO. 

SYlrO. 

SX2=0. 

SX3=0. 

SX4=0. 

SXIYUO. 

SX2Y1=0. 

XN=2*NFIT+1 

C 

C  COMPUTE  TIE  MATRIX  ELEMENTS 
C 

DO  20  M=-NFIT,NFIT 

IF  (N+M-l.LT.l)  GO  TO  25 
IF  (N+M-I.GT. NDATA)  GO  TO  25 
YUY(N4M-1) 

25  XUFLOAT(m) 

X2=X1»X1 

X3=X2»XI 

XA=X3»X1 

X1YUX1»Y1 

X2Y1=X2«Y1 

SYlrSYUYl 

5XUSXUX1 

SX2=5X2+X2 

5X3=SX3*X3 

SXA=SXAtXA 

SXlYlzSXlYltXlYl 

SX2YUSX2YUX2Y1 

20  CONTINUE 

MATRIX(l,l)rXN 
MATRIX(2, 1)=5X1 
MATRIX(3,1)=SX2 
MATRIX(1,2)=SX1 
MATRIX(2,2)=SX2 
MATRIX(3,2)=SX3 
MATRIX(1,3)=SX2 
MATRIX(2,3)=SX3 
MATRIX(3,3)=SXA 
MATRIX(1,4)=SY1 
MATRIX(2,A)=SX1YI 
MATRIX(3,A)=SX2Y1 


c 

C  CHECK  ^^€  SIZE  OF  THE  MATRIX  ELEMENTS 
C 

DO  30  1=1,2 

DO  AO  J=I,3 
11=1 
J1=J 
C 

C  IF  THE  MULTIPLIER  WOULD  BE  >  1,  PIVOT  Tit  MATRIX 
C 

IF  (ABSIMATRIX  (I, I) ) .LT . ABS(MATRIX( J, I ) 
1  ))CALL  INVERT(I1,J1, MATRIX) 

AO  CONTINUE 

C 

C  COMPUTE  THE  REDUCED  MATRIX 
C 

DO  50  J=l+1,3 

MULT=-MATR1X( J, I )/MATRIX( 1,1) 

DO  60  K=1,A 

MATRIX! J,K)=MATR1X(J,K)+ 

1  MULT*MATRIX{I,K) 

60  CONTINUE 

50  CONTINUE 

30  CONTINUE 

C 

C  CALCULATE  THE  VALUES  OF  THE  COEFFICIENTS 
C 

C(N)=MATRIX(3,A)/MATRIX(3,3) 

B(N)=(MATRIX(2,A)-C(N)*MATRIX(2,3))/MATRIX(2,2) 

A(N)=(MATRtX(l,A)-B(N)«MATRIX(l,2)-C(N)»MATRIX(l,3)) 

1  /MATRIX! 1,1) 

10  CONTINUE 


RETURN 

END 


SUBROUTINE  INVERT! I, J,MATRIX) 

»•••••»•«•••«•••»•«»»«>«•«•««»«  «  >«•»•«••»«»*••»•*««••» 

C 

C  THIS  SUBROUTINE  PIVOTS  THE  MATRIX  TO  REDUCE  ERROR  PROPAGATION. 

C 

DIMENSION  MATRIX!3,A) 

DOUBLE  PRECISION  MATRIX,  SMALL 
DO  10  K=1,A 

SMALL  =  MATRIX!l,K) 

MATRIX!I,K)  =  MATRIX! J,K) 

MATRIX!J,K)  =  SMALL 
10  CONTINUE 

RETURN 
END 


0*#*****#»#»*##*»»**»***#*tt#***#lunni#*»##«  ««»#*»*«»*««#»#»««*«*»##*■******«*** 

SUfjROUTINE  FILTER!Y,A,B,C) 

C 

C  THIS  SUBROUTINE  MAKES  A  PREDICTION  FOR  TIE  NEXT  POINT  BASED  ON  THE 
C  COEFFICIENTS  CALCULATED  FOR  TIC  FOUR  POINTS  AROUND  IT. 

C 

DIMENSION  Y!NDATA),A!N0ATA),B!NDATA),C!NDATA) 

COMMON  NOATA,  NFIT 
DOUBLE  PRECISION  Y,A,B,C 
IEND=NDATA-2 
C 

C  THE  END  POINTS  HAVE  A  DIFFERENT  WEIGHTING 
C 


Y(l)=.6»A!l)t.3»!A!2)+B!2)tC!2))+.l»!A!3)t2*D!3)tA«C!3)) 
Y!2)=.2*!A!l)-B!l)tC!l))+.5»A!2)t.2»!A!3)+B!3)4C!3))4.1«!A!A)+2 
1  B!A)4A»C!A)) 


WRITE  (21,100)  Y(I) 

10  CONTINUE 

CLOSE  (21) 

NBUEF=NBUFF+1 
RETURN 

100  FORMAT  (020.12) 

END 

Q********************************«****»*»*»**«************«**»«, 

SUBROUTINE  RBUFFl(Y) 

««««*»»««« «»«•«*««»•««•«• 

c 

C  THIS  SUBROUTINE  ADDS  AN  EXPONENTIAL  DECAY  TO  T^€  BEGINNING  OF  THE  SPECTRA 
C  SO  THAT  THE  FILTERED  SPECTRA  DO  NOT  GET  A  BIG  IMPULSE  AFFECT  ON  Th€  ENDS 
C 

DIMENSION  Y(120) 

COMMON  NDATA.NFIT.NBUFF 
DOUBLE  PRECISION  Y 
DO  10  1=110,1,-1 

Y(I+10)=Y(I) 

10  CONTINUE 

C 

C  COMPUTE  THE  FIRST  POINT  BASED  ON  THE  SLOPE  AND  END  POINT  IN  THE  SPECTRUM 
C 

FIRST=2.*Y(11)-Y(I2) 

DO  20  1=10, 1,-1 

XPOW=(-(10.-I)*(10.-1)/50.) 

Y(I)=FIRST*EXP(XPOW) 

20  CONTINUE 

RETURN 
END 

»»»*»«»*««*»»«»**«»«*«»»«»»«»#»««••««•«»««««««*««««*««««««*««««  ««««*«««•»« 

SUBROUTINE  RBUFF2  (Y.INT) 

c 

C  THIS  SUBROUTINE  BUFFERS  THE  LAST  DATA  SET 

I  C 

I  DIMENSION  Y(120) 

COMMON  NDATA,NFIT,NBUFF,NFLAG 
DOUBLEPRECISION  Y 
NFLAG=1 
C 

C  THE  LAST  POINT  IS  BASED  ON  THE  SLOPE  AND  END  POINT  OF  THE  DATA  SET 
C 

YLAST=2*Y(110)-Y(109) 

DO  10  1=1,10 

XP0W=-I»I/50. 

Y(110+I)=YLAST*EXP(XP0W) 

10  continue 

RETURN 
END 

SUBROUTINE  ERROR(IOS) 

IT  «#«##««**#»»»#»«»»»»•»»»»»««»»»»»»#»»#««»»«»»«  •««««»»»»»««»»«»« 

C 

C  THIS  SUBROUTINE  IS  PURELY  PRECAUTIONARY,  IT  WARNS  IF  AN  I/O  ERROR  IS 
C  ENCOUNTERED. 

C 

WRITE  (5,100)  lOS 

100  FORMATC  I/O  ERROR  ENCOUNTERED.  I/O  ERROR  IS  M3) 

RETURN 
END 

»»»««•*•»«»»•*««»•»•«•«»  »«••••••«•• 

SUBROUTINE  MINUS(riLEl ,FILE2, INT ) 

»»««*«» ft »»«•«*»»«««•»••«•«»»»« »««««« 

c 

C  THIS  SUBROUTINE  WILL  READ  THE  TEMPORARY  FILES  AND  TI€N  SUBTRACT  THEM 
C  POINT-BY-POINT 
C 


DIMENSION  Y1(IOOO),Y2(IOOO),FILE1(6),FILE2(/;),STOR1(2),STOR?(?) 


Y(1EN0+1)=.1*(A(IEND-1)-2*B(1EN0-1)+4*C(IEND-1))  +  .2*(A(1EN0)- 
1  B(IEND)+C(IEND))+.5*A(IENO+l)+.2*(A(IEND+2)+B(IEND+2)+C(IENI>2)) 

Y(IEND+2)=.6*A(lEND+2)  ♦  .3*(A( IEN0+1)-B( IEND+1)+C( IEND+1) )  + 

1  ,l*(A(IEN0)-2*B(IEND)t4»C(lEND)) 

C 

C  MOST  POINTS  ARE  WEIGHTED  I€RE  1,2, 4, 2,1 
C 

DO  10  I=3,IEND 

Yl  =  .l»(A(I-2)-2.»B(I-2)+4.<'C(I-2)) 

Y2=.2*(A(I-1)-B(I-1)+C(I-1)) 

Y3=.4*(A(I)) 

Y4=.2*(A(I+1)+B(I+1)+C(I+1)) 

Y5=.1»(A( I+2)+2.*B( 1+2)44. »C(  1+2)) 

Y(I)=Y1+Y2+Y34Y4+Y5 
10  CONTINUE 

20  FORMAT  (020.12) 

RETURN 

END 

«««•«««»«««««»««««  »4  »*»««»«*«»«**» 

SUBROUTINE  RBUFF(Y,FILENA,INT) 

■  «««««*«#«  c  ))««««««*««««««««««««««« 

c 

C  THIS  SUBROUTINE  READS  A  100  POINT  PORTION  OF  THE  SPECTRUM  AND  BUFFERS 
C  IT  BY  10  POINTS  ON  EITHER  SIDE. 

C 

DIMENSION  Y(120),FILENA(6) 

COMMON  NDATA,NFIT,NBUFF,NFLAG 

DOUBLE  PRECISION  Y 

INT=0 

OPEN  (20, NAME=FILENA,ACCESS= ’DIRECT ’.RECL =22) 

11=NBUFF *100-9 
IEND=I1+119 
DO  10  1=11, lEND 

IF  (I.LE.O)  GO  TO  10 

READ  (20,100,REC=I,ERR=2O,IOSTAT=I0S)  Y(INT+1) 

INT=INT+1 

10  CONTINUE 

C 

C  CREATE  A  LEADING  BUFFER  IF  THIS  IS  THE  FIRST  DATA  SET  IN  THE  SPECTRUM 
C 

IF  (NDUFF.EQ.O)  CALL  RBUFFl(Y) 

CLOSE  (20) 

RETURN 

20  IF  (lOS.NE.SlO)  CALL  ERROR(IOS) 

C 

C  CREATE  A  TRAILING  BUFFER  IF  THIS  IS  THE  LAST  DATA  SET  IN  THE  SPECTRUM 
C 

CALL  RBUFF2(Y,INT) 

CLOSE  (20) 

100  FORMAT  (D20.12) 

RETURN 

END 

*»»*»»#»»  »»#»»*«»#«»«»»»  ■»*)HHH»»4«********* 

SUBROUTINE  WBUFF  ( Y.FILENA, INT ) 

»»•»*»* ►*»**»»»»**»»*••»*•»»••••»«••»•»»»•*««»»♦»• *«»»»»»»»»***»**•* 

C 

C  THIS  SUBROUTINE  TRUNCATES  THE  OUTER  10  POINTS  ON  EACH  SIDE  AND  WRITES 
C  Tie  REMAINDER  TO  DISK  IN  SEQUENTIAL  ORDER. 

C 

DIMENSION  Y(120),ST0RE(2),FILENA(6) 

COMMON  NDATA,NFIT,NBUFF,NFLAG 
DOUBLE  PRECISION  Y 
ST0RE(1)=FILENA(2) 

ST0RE(2)='.TMP' 

Il=N0in  •100*11 
IEND=Il+99 
J  =  ll 

OPEN  ( 2 1 , N AME =S TORE , ACCE SS= ' APPEND ' ) 

no  10  1=11,110 


DOUBLE  PRECISION  Y1,Y2 
COMMON  NDATA,NFIT,NBUFF,NFLAG 
ST0R1(1)=FILE1(2) 

STORI(2)='.TMP’ 

ST0R2(1)=FILE2(2) 

ST0R2(2)='.TMP' 

OPEN  (20,NAME=ST0R1 ,ACCESS= 'SEQUENTIAL ' ) 
OPEN  (21 ,NAN€=ST0R2 , ACCESS= 'SEQUENTIAL ' ) 
READ  (20,100,END=20)  Y1 
READ  (2I,100,ENO=30)  Y2 
IEND=NBUFF*100 
REWIND  (20) 

DO  10  I=1,IEND 

Y1(I)=Y1(I)-Y2(I) 

CONTINUE 
WRITE  (20,100)  Y1 
CLOSE  (20) 

CLOSE  (21) 

RETURN 

FORMAT  (D20.12) 

END 


^  <.*r 
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